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ABSTRACT 

Formulating the equations of motion for cosmological bodies (such as 
galaxies) in an integral, rather than differential, form has several advantages. 
Using an integral the mathematical instability at early times is avoided and 
the boundary conditions of the integral correspond closely with available data. 
Here it is shown that such a least-action calculation for a number of bodies 
interacting by gravity has a finite number of solutions, possibly only one. 
Characteristics of the different possible solutions are explored. The results are 
extended to cover the motion of a continuous fluid. A method to generalize an 
action to use velocities, instead of positions, in boundary conditions, is given, 
which reduces in particular cases to those given by Giavalisco et al. (1993) 
and Schmoldt & Saha (1998) . The velocity boundary condition is shown to 
have no effect on the number of solutions. 



Subject headings: cosmology:theory — galaxies: kinematics and dynamics- 
galaxies: formation — methods: numerical 
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1. Introduction 

The present motions of cosmological objects, in particular galaxies, are functions of 
their past history. In principle one might discover the shape of the past by calculating 
presently observed positions and motions backward. However, in doing this we are 
faced immediately with two problems. First, their velocities in the plane of the sky are 
not known, and their distances not known accurately; so perhaps half the information 
needed to start the calculation by Newton's equations is there. Second, if the trajectories 
of the galaxies are to be traced back to very early times distances become very small 
and corresponding gravitational forces very large. Small errors in present velocities or 
positions become heavily magnified, resulting in galaxies being formed at infinite speeds. 
The problem is mathematically unstable, rather like trying to roll a marble to the top of 
a glass mountain, and requiring that it stop exactly on the summit^. 

To avoid these difficulties Peebles (1989, 1990, 1994) formulated the problem in 
integral rather than differential form. This traded the relative simplicity and definiteness 
of differential equations for the stability of the integral. The most important consideration 
in moving from the differential to integral form of the problem (apart from the mechanics 
of implementation) is the fact that, with the same boundary conditions, an integral 
calculation may produce several (or many) solutions. An obvious question to answer is 
just how many there are. This is something more than a purely mathematical concern. 
Of course, if the numerical calculation of solutions can be guided in some way there is 
the potential for a large savings in computer time, and if the number of solutions is 
limited the search may be stopped when all are found. Conversely, if the number of 
solutions is very large or infinite, the usefulness of the calculation is thrown into doubt 
(unless some method of selecting more probable solutions is found). But the question is 
more fundamental than that, for the variational formulation of the cosmological problem 
corresponds closely to the limits of our knowledge. When the present radial velocities and 
positions on the sky of a number of bodies are specified and the Big Bang postulated, 
we find the end conditions are fixed; the action is determined by relevant physics. The 
mathematical question is thus transformed into a cosmological one. 

The subject of this study is the mathematical theory of variational calculations 
as applied to the cosmological problem. That problem is defined as the determination 
of the motion of a number of bodies moving under gravitational interaction, with the 
requirement that all bodies must be at the same point (in proper coordinates) at t — 0. 



1 Valtonen et al. (1993) have found some possible solutions for the motion of the major 
galaxies in the Local Group and the Maffei 1/IC 342 Group by integrating equations of 
motion forward from an early time. However, it is not clear that this method is generally 
applicable, and in any case requires a great deal of hunting about in parameter space; for 
their result, the Valtonen group integrated ten thousand situations. 
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Newtonian, rather than relativistic, calculations are employed throughoutf]. 

The cosmological problem may be interpreted as a rough approximation of the motion 
of galaxies, each galaxy simulated by a point mass interacting only through gravity. This 
is the way most least-action calculations have proceeded, and is not a bad approximation 
considering the uncertainties in such data as distances and masses. It would be more 
accurate, however, to consider the objects to represent the dark matter halos of galaxies 
(which as far as is known interact only through gravity). The point-mass approximation 
provides a reasonable simulation of gravitational effects, since multipole moments decay 
rapidly with distance (Dunn & Laflamme 1995 found them to be quite unimportant), 
and in any case the conclusions of this study are not affected by the detailed form of the 
gravitational potential used. 

Of course, identifying whole galaxies with single bodies ignores internal structure 
(which may be significant in some cases) and the effects of mergers (which certainly 
are significant); Dunn & Laflamme (1995) found some additional problems. To address 
these matters one must turn to a continuous fluid formulation of the problem. Section 6 
generalizes the discrete-body results to this more complicated situation. 



2. The First Variation 



Consider first the problem of minimizing the integral 



/= [ tl (T-V)dt= f 1 L( qi ,qi,t) 
Jt Jt 



dt 



(1) 



where the kinetic energy T is quadratic in the generalized velocities q\ (or, alternatively, 
in the generalized momenta dL/dq\ = and the potential V does not depend on 
velocity. The end points qi(t ) and qi(ti) are given. (This is interpreted dynamically by 
constructing a path r(t) in 3n-dimensional space using the vectors rj(qi(t)), j = 1 to n, 
i = 1 to 3n, where the Vj are the paths of the n bodies in 3-dimensional space.) For small 
variations the change in the action is given by a truncated expansion in a Taylor's series 
(treating qi,(ji as independent variables): 



8 / L{q il q h t)dt = / }2 

Jto Jto • 



dL r dL r . 
oqi oqi 



dt = 0. 



(2) 



Equation (0) can be integrated by parts to give 
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dqi 
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dqi 



(3) 



to 



2 See Peebles (1980) and Bondi (1960) for the validity of this approach. 
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and if the variation in r vanishes at the end points (that is, if the end points are fixed) the 
boundary term is zero. The requirement that the integral vanish for arbitrary variations 
Sr results in the Euler-Lagrange equations: 



dL d (dV 
dqi dt \dq i/ 



0. 



(4) 



These are the dynamic equations, identical with those derived from (for example) forces 
and accelerations. The correspondence between the dynamic equations and the vanishing 
of the first variation of equation (Jl|) is Hamilton's Principle. 

A path which minimizes the integral will thus always satisfy the dynamic equations. 
However, the converse is not necessarily true: a path satisfying the Euler-Lagrange 
equations is not guaranteed to provide a minimum of the corresponding integral. For 
sufficiently small path lengths a minimum does result (see, for example, Whittaker 1959, 
pp. 250-2); beyond a certain point the path will make the integral stationary, but not 
necessarily a minimum. Finding the point that determines the limit of application of the 
least-action technique (strictly interpreted) to the dynamical problem will be discussed 
belowf]. 

Adding a total derivative (in multiple dimensions, a divergence expression) will not 
change the Euler-Lagrange equations (see Courant and Hilbert 1953, p. 296) but can 
change the boundary terms. For instance, varying 



t 



L 
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0. 



(5) 



(6) 



If the original Lagrangian is quadratic in <jj, recovery of the Euler-Lagrange equations 
requires that either qj = or dqj = in each boundary term. In the second case, it 
is the velocity at the end point which is the fixed boundary condition, rather than the 
position. This raises the possibility of using a radial velocity, rather than a distance, 
as the end point in a cosmological calculation. In fact Giavalisco et al. (1993) have 
considered such a mixed boundary condition, in one place using it to modify a set of 
approximating functions and in another expressing it as a canonical tranformation of 
variables. Their approaches are in practice equivalent to this one. Schmoldt & Saha 
(1998) have succeeded in using a velocity endpoint in their numerical calculation. 



It is straightforward to show that the boundary term added here does not change 
Whittaker's conclusion above. Note that the coordinates qj which provide velocity 



3 Whittaker calls Euler's Principle, which appears later, the Least Action Principle; 
however, Hamilton's Principle has also been called this. To distinguish between the two I 
will use the names of the mathematicians and call them collectively Least Action Principles. 
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boundary conditions may be all or only some of the total number of coordinates g^, as 
long as the total derivative is adjusted accordingly 



2.1. Variable Endpoints 

If the radial velocity is to be used as a form of endpoint, rather than the (less 
accurately known) radial distance, the theory of "free" endpoints (constrained to move 
on a manifold of some description) comes into play. In addition to the Euler-Lagrange 
equations, the solution must now satisfy the transversality condition which results from 
the minimization of the variation at a free end pointQ. As expressed in Morse's (1934) 
notation, this condition is 

(.-pf^ + Ei^O (7) 

where the superscript s denotes a differential taken along the end manifold (s takes on 
values designating the initial or final end points) and L is the integrand. 

Applied to the Lagrangian for a number of bodies moving under their mutual gravity 

L = J2 m i ( ~ (jl + r 2 M + A sin 2 0$) + G £ (8) 
i \ z j<i \ r v\J 

and fixing the time, the transversality condition is 

]T rrn (fidr? + rjOid&t + r\ sin 2 difadtf) = (9) 

i 

or more compactly 

P • dr s = (10) 

where P is the total momentum and dr s any vector in the end manifold, a general result 
constraining the end manifold. If the velocity action, expression (||), is used a modified 
form of the transversality condition applies: 

mi (fidrf + rfdidBl + r] sin 2 Oifad^f) - £ m^idr* - £ m^df* = 0. (11) 

i i i 

If the end point under consideration has fixed angles and radial velocities, the left hand 
side of equation ( |TT] ) vanishes identically. The velocity-action transversality condition 
thus tells us nothing about the manifolds on which the end-points lie. At the same time, 
the velocity action imposes no additional restrictions on the end manifolds over the 
position action. 



4 See Appendix A for detailed formulae. 
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3. The Second Variation 

We now come to the question of how far the minimization of the action integral can 
be used to reproduce the dynamic equations, that is, the limit of the least-action method 
strictly defined. The limit may be pictured geometrically by using kinetic foci as defined 
by Thompson and Tait (1896, section 357, p. 428)0: "If, from any one configuration, two 
courses differing infinitely little from one another have again a configuration in common, 
this second configuration will be called a kinetic focus relatively to the first: or (because 
of the reversibility of the motion) these two configurations will be called conjugate kinetic 
foci." It can be shown (for instance, by Whittaker 1959, pp. 251-3) that the action is 
neither a maximum nor a minimum over a path which includes a pair of kinetic foci. 

More intuitively, if two paths infinitesimally close to each other between the same 
pair of end points both satisfy the Euler-Lagrange equations, the action (along either 
path) can no longer be a minimum (and the variation of the variation between them must 
vanish) . 

It is easiest to picture kinetic foci using a toy dynamical problem, that of finding the 
motion of a ball rolling on a large sphere, without friction or other complicating effects. 
Geodesies on the spherical surface are paths of least action in this case. Clearly there is 
a unique minimum path for points close together; that is, if two end points are chosen 
near each other, a portion of the great circle joining them gives the least action. As the 
points are taken farther and farther apart the action increases while staying a minimum. 
When the points are taken to be opposite each other, however, there is an infinite number 
of solutions all of the same length. Dynamically, the ball leaves the starting point and 
follows a geodesic to the antipode; but if it had left the starting point at a slightly 
different angle, it would still pass through the same antipode. On a sphere, then, kinetic 
foci are exactly opposite each other. 

Still considering the motion of the ball dynamically, if the trajectory is extended, 
worse happens. The path taken, which passes more than halfway around the sphere, is 
actually longer than paths which follow the complement of the great circle. In fact it is 
longer than some small circles. 

Geodesies on a torus provide greater complexity. Given any two points on the torus 
there will be a global minimum; a local minimum path, going the other way around the 
major radius; and an infinite series of other local minima, wrapping around one leg or the 
other of the torus. None of these can be continuously varied into another because of the 
different number of wrappings. The action (the length of the geodesic) tends to infinity 
as the wrappings increase. 



5 Page and section numbers are identical in the 1962 Dover reprint. 
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A more mathematically rigorous and useful, but also more complicated, treatment of 
the second variation takes us into Morse Theory. 



3.1. Morse Theory 

Marston Morse (1934) conducted an extensive study of the general topological 
properties of variational problems and their solutions. A summary of some of his results 
is presented belowf]. 

An extremal is a path which satisfies the Euler-Lagrange equations. A critical 
extremal is one which makes the action a minimum. 

A function which will play a part in what follows is defined by 



d 2 L 2 n d 2 L . <9 2 L. ? 

™ M= M s < +2 ^ iSi+ m s * (12) 



for some integrand L and functions Si(t). The characteristic form is[] 

Q(z, A) = f 1 2 {Q( Si , - Xs lSj ) alt. (13) 

J to 

For a given extremal with fixed end points, the accessory boundary problem is defined as 

d fdfl\ dQ , 

*U)-^ + Asi = a (14) 

A solution Si to this equation not identically zero is an eigensolution (sometimes 
eigenvector) and A an eigenvalue. The index of an eigenvalue is the number of linearly 
independent eigensolutions corresponding to the eigenvalue. 



For free end points the characteristic form is 

rti 

bhkUhUk + I 

h,k 

where bhk, Uh and Uk are derived from the second variation at the end points and are 



rti 

Q(z,X) =2Zb hk u h u k + I 2 (£l{si, Si) - XsiSj) dt (15) 

h h J to 



given in Appendix A. The accessory boundary problem, equation (0), stays the same in 



form but the solution s, must now satisfy the transversality condition. 



6 A shorter and somewhat more accessible presentation of most of Morse's results is 
found in Milnor (1963). 

7 Here z is used as a shorthand symbol for the collection of functions Sj, and below it 
includes also Uh and u k . 
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If the problem is to find the geodesic in a space of a given metric between two 
manifolds of some description, the coefficients bhk are a measure of the curvature of the 
manifolds. In particular, when the coefficients vanish the manifolds are flat. 

If A = the accessory boundary problem becomes the Jacobi equation (not to be 
confused with the Hamilton- Jacobi equation), which is identical to the perturbed Euler- 
Lagrange equation. If there exist two points on an extremal at which an eigensolution with 
eigenvalue zero vanishes, these points are conjugate points. A non- degenerate extremal is 
one which has no zero eigenvalues in the accessory boundary problem. It is easily shown 
that conjugate points (mathematically defined) and kinetic foci (dynamically defined) are 
identical. 

Determining the least-action limit for a given variational problem is thus the same 
as determining the first zero of the Jacobi function (after the initial point). This 
determination is not generally an easy thing to do. To take a specific example, for a group 
of bodies moving under each other's gravity the Jacobi function s satisfies 



Clearly, solving this is not a convenient way of finding kinetic foci. Not only is this less 
amenable to integration than the original dynamic equation, but the original equation 
must be solved first (which makes the locating of kinetic foci as a step in solving the 
dynamical problem rather pointless). A more practical method for use in calculation is 
called for. 



Choquard (1955) studied the motion of bodies in strongly anharmonic potentials in 
the context of a semi-classical treatment of Feynman integrals. He found that multiple 
solutions to a dynamical problem were possible through the action of "forces of reflection", 
which allowed indirect paths from one end point to the other. In an indirect path, which 
corresponds to a stationary rather than a minimum action, at some time between the end 
points 




(16) 



3.2. Choquard's Criterion 



d_ 

dt 



(T) = 




(17) 



That is, the momentum must be normal to the force. 



To make this reasoning directly applicable to the problem at hand, consider a solution 
to the dynamical equations with a given set of end points, r(£); it must conserve total 



-10- 



energy E, made up of a kinetic part T and a potential part V. A varied path r(i) + s(t) 
(where s(t) is a Jacobi function) also conserves an energy E + 5E = T + 5T + V + SV, 
and thus the Jacobi function itself conserves ST + SV. Since V is a function only of r, not 
r, 5V is a function only of s, and for |s| small (which it is by definition) a linear function. 
This means that 5V reaches its extreme value when s does, and at the same point ST has 
an extremum. Since |s| can vanish only after its maximum, this point of extremum must 
occur before a conjugate point. The extremum of ST is given by 

±(ST )= 0. (18) 

Writing kinetic energy in terms of momentum, 

T + ST = -L( p + 5p) 2 
2m 

ST ~ — p-5p 



m 

= — p-(-WUt (19) 
m 

so the condition for an extremum of the variation in kinetic energy is 

= 

-p-(-vy) = o (20) 

m 

and Choquard's criterion is recovered. For a situation with multiple particles, the varied 
path s may be taken to be different from zero for only one of the particles. This leads to 
the conclusion that a conjugate point may occur only after the point where the momentum 
is normal to the force on some body in the system. Keeping in mind the identity of 
conjugate points and kinetic foci as well as Whittaker's result (above), Choquard's 
criterion gives a lower bound to the applicability of the least-action calculation. Following 
the trajectory of a dynamic system from the initial point, it is a minimum of the action 
at least until the momentum of some body is normal to the force on that body. This 
provides some insight into the shape of stationary-solution trajectories, as well as (with a 
further result of Morse, below) allowing a conclusion to be drawn as to the total number 
of solutions of all kinds^. 



3.3. More Morse 



There are several more results from Morse (1934) which are of use in the present 
problem. First we require a few more definitions: 



8 Choquard (1955) notes that his criterion does not apply to situations in which 
the trajectory is always normal to the acceleration, as in (for example) circular motion. 
However, these situations are generally symmetrical enough to allow the useful application 
of Jacobi functions. 
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A Riemannian space possesses a positive-definite metric which can be expressed as a 
quadratic form: 

ds 2 = ^2g ij dx l dx j . (21) 

i,j 

The connectivity Pk of a spacef] is the number of distinct homologous families of 
figures of dimension k + 1; that is, within each family one figure can be transformed 
into another by a continuous transformation, but a figure in one family cannot be so 
transformed into a figure in another. On a sphere, for instance, the connectivity Pq is 
one, since any line may be transformed into another by a continuous transformation. On 
a torus Pq is infinite, since there is an infinite number of families of curves distinguished 
from each other by the number of times they wrap about the large or small radii. 

Morse is concerned with the connectivities of the functional domain Q of admissible 
curves for a given variational problem^, that is those curves which have the required 
end points and are continuous along with their first derivatives. For the case of a set 
of trajectories in three dimensional space it is easiest to consider them transformed 
into a single trajectory in 3n-dimensional space, between two end points representing 
the starting and ending configurations. Each point in Q represents a trajectory in the 
3n-dimensional space. Since no points in 3n-space are excluded, any trajectory can 
be continuously transformed into any other; so any point in Q can be continuously 
transformed into any other. Any line in Q can then be transformed point by point into 
any other line, any plane figure likewise, and so on for all dimensions. Consequently each 
connectivity of the space of trajectories is one. 

Morse's important results are: 

An extremal which affords a minimum has no negative eigenvalues in the associated 
boundary problem. This is equivalent to saying it contains no conjugate points. Further, 
the number of conjugate points of an end point of an extremal g on g is equal to the 
number of negative eigenvalues in the associated boundary problem. 

The index of an extremal is the sum of the indices of the conjugate points of an end 
point on the extremal. 

The conjugate points of an end point of an extremal g on g form a set of measure 
zero. This means they are isolated (and thus much easier to deal with). More importantly, 
it means that the probability of choosing a pair of conjugate points by chance when 
setting up the variational problem is essentially zero. 



9 This is not to be confused with the connection of a space, or whether a space is simply 
connected (themselves distinct topological concepts). 

10 This is not the function Q(s, s) found above and in Appendix A. The ambiguity in 
notation is regretted, but it should not lead to confusion. 
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If for a given Riemannian space R and terminal manifold Z there exists an integral 
I defined on R such that all critical extremals are non-degenerate, then the number 
of distinct extremals of index k is greater than or equal to the connectivity Pk of the 
functional domain Q. If the extremals are of increasing type, the number of extremals of 
index k is equal to the connectivity P}~. 

The last is a most useful result. However, to apply it we must show that the 
variational problem meets the requirements. 

As demonstrated for example by Whittaker (1959, pp. 247-8, 254)0 the dynamic 
equations of a system which has an integral of energy E can be derived by requiring that 
the variation of the integral 

r 2Tdt (22) 



(where T is the kinetic energy) vanish, for a fixed value of E. This formulation is known 
as Euler's Principle. For a system in which the total energy E is the sum of the kinetic 



energy T (quadratic in velocities) and potential energy V, the integral (|22| ) can also be 
written as 



J 2(E-V) 1/2 (T) 1/2 dt 

2 (E - V) 1/2 (ayiV') 1/2 dt. (23) 



This is the integral giving arc length on a surface of metric 

g ij = 4(E-V)a ij (24) 

and the space of the possible solutions to our variational problem when formulated this 
way is indeed Riemannian, with geodesies corresponding to variational solutions. The 
cosmological variational problem in proper coordinates can be put into this form, so 
Morse's results apply. Further, the integral is of increasing type (since kinetic energy is 
always positive), so the number of solutions is given definitely. 



Unfortunately, Morse's result cannot be applied directly to Hamilton's principle, as 
the functional domain is not Riemannian. The Lagrangian cannot be put in the necessary 
form of expression fl23|), since the potential energy forms a separate term which is not 
quadratic in the velocity differentials. This means that, if the preceding theory is to be 
used, either the problem must be put in Euler's form or some connection of a topological 
nature must be made between the two Least Action principles. The latter is addressed in 
the next section. 



Before leaving Morse Theory, however, it is worthwhile to note the effect of the end 
manifolds of the velocity action on the number of solutions. Comparing the characteristic 
forms and accessory boundary problems between fixed-point and manifold situations, we 



n See also Arnold (1989), pp. 245ff. 
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find that the difference lies in the transversality condition and the quantities bhkUhUk- It 
has already been shown that, using Euler's Principle, the transversality condition holds 
identically. 

For an initial configuration manifold and Hamilton's Principle, combining two 
formulae from Appendix A, 



hkUhUk 



n ^(dLdt s dqt dLd 2 q^ 2 




y dqi de de dq\ de 2 / { 
= ^ m \^ + r 2 e i -^ + r 2 ^ 2 e i ^). (25) 

Comparing this with equations (||) and ( |10D the manifold curvature expression bhkUhUk is 
seen to be the dot product of the total momentum vector with a vector in the end-manifold 
surface: 

d 2 

bhkU h Uk = P • ^jx s (26) 

which is, near an extremal, 

bhkUhUk = P • Ax s /e 2 . (27) 



For an extremal, the tansversality condition requires that P • Ax s = (see equation |10|) . 
More generally, since P is a constant of motion belonging to the solution (not related to 
any particular variation around it) and e and Ax s are independently arbitrary, P • Ax s 
cannot vary with e, and thus must vanish (unless bhkUhUk is allowed to be infinite, a 
pathological case I propose to ignore). The dot product vanishing causes bhkUhUk to 
vanish as well. The same result holds if Euler's Principle is used. 

For the velocity action and the final end manifold, with angles and radial velocities 
fixed, bhkUhUk vanishes identically using either Least Action Principle. Thus for each case 
the fact that the ends of the action integral lie on manifolds and not on fixed points is 
irrelevant to the number of solutions. Interpreted geometrically, the manifolds are flat 
surfaces. 



4. The Number of Solutions 

4.1. Catastrophe Theory 

Applying Morse's result on the number of solutions to the problem formulated using 
Euler's Principle (which is the only way it is directly applicable), we find that there is one 
minimal extremal, plus one non- minimal extremal for each integral number of conjugate 
points. There are values of total energy for which there are no solutions. Most obvious 



-14- 



are those below the potential energy of the final end point; those are excluded from 
the functional domain at the outset. For values of the total energy in a gravitational 
system which are positive, especially strongly so, there can be only one solution since the 
Jacobi function for a nearly straight-line trajectory never returns to zero; the saddle-point 
solutions for these situations can be thought of as occurring at infinite values of total 
time. 

But the calculation contemplated (and as performed by Peebles and those following 
his technique) uses Hamilton's Principle. To apply Morse to Hamilton a connection 
must be made between them. This can be done by way of the dynamical equations and 
Catastophe Theory. 

Consider a two-dimensional slice of the functional domain Q, the dimensions being 
the Eulerian action (time integral of the kinetic energy T) and the Hamiltonian action 
(time integral of the Lagrangian function T — V). Choose the slice so that it contains all 
the extremals of the problem (see Figure |IJ). A given value of total energy will plot as a 
curve in this slice, with a minimum of / Tdt at the location of the least-action trajectory 
and other extremals spaced along it (the latter may show as maxima, minima or points 
of inflection in this plot). Since the integrand of Euler's Principle is positive-definite, the 
least-action solution takes the mimimum time of all the solutions for a given energy; the 
saddle-point solutions take increasing amounts of time for increasing index. All solutions 
will be equilibrium points for the potential represented by the action. The slice is thus a 
Poincare diagram to which Catastrophe Theory applies, with total energy as the control 
parameter. The minimum solutions correspond to stable equilibria, the saddle-points to 
unstable equilibria. 

In the same slice plot curves of constant total time. The index of a given extremal 
depends only on the number of kinetic foci of the trajectory, not on the action principle 
(if any) used to calculate it; in addition, all extremals are solutions to the dynamic 
equations. Thus the minimum of / Ldt for each time corresponds to a minimum of 
/ Tdt for fixed total energy, and the non-minimum extremals will similarly correspond to 
non-minimum extremals of the Eulerian action for other values of total energy. Again, we 
have constructed a Poincare diagram (rotated 90° with respect to the first), with total 
time as the control parameter. 

There is at least one least-action solution for a given value of time. If there were 
two (or more), the chain of Eulerian least-action solutions would have a maximum or a 
minimum in total time, as shown in Figure |l|. There are several reasons why this cannot 
happen; two are outlined below. 

First, viewed as a Poincare diagram in Hamiltonian extremals, Figure [I] requires 
two chains of similar (stable or unstable) equilibria to meet. This sort of topology, 
a bifurcation without an exchange of stability, is forbidden by Catastrophe Theory; 
therefore there is only one least-action Hamiltonian extremal. Similarly, there can only 
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be one saddle-point Hamiltonian extremal for each saddle-point Eulerian extremalQ 
Second, note that at a bifurcation point (point C in figure [l|) 



for the action / and variables q in the two-dimensional slice; but this is just the 
requirement for a degenerate extremal, to which Morse's results specifically do not apply. 
Since, as noted above, these require some special symmetry in the problem and are almost 
impossible to generate by chance, it is reasonable to assume that our problem does not 
have themj^|. 

We are finally in a position to determine how many solutions there are to the 
cosmological variational problem. 

If the question is posed in a strictly proper-coordinate, Newtonian manner it comes 
out something like this: given a number of bodies moving under the influence of each 
other's gravity, all constrained to occupy the same position at time zero, and having given 
positions (or positions in two dimensions, radial velocities in the third) now; how many 
possible trajectories are there? 

If the problem is formulated using the Eulerian action (minimum kinetic energy 
for fixed total energy), the space of solutions is Riemannian and the extremals are of 
increasing type. There is thus one minimum (and one stationary solution for each number 
of kinetic foci). By way of Catastrophe Theory this is connected to the Hamiltonian 
action (the form in which the question is asked above), which excludes some solutions 
which require a different value of total time. There is one minimum solution and a finite 
number of stationary solutions. 

For small values of total time the energy will be forced to be positive (in order for 
the system to get from one configuration to the other, the speeds must be large, hence 
the kinetic energy large and positive) and only the least action solution will appear. This 
idea will be expanded below. 

Note that if there is no integral of energy these results do not apply. Thus 
if a calculation attempts to compute the trajectories of a number of galaxies in 



12 Expositions of Catastrophe Theory are found in Lamb (1932, sect. 377, pp. 710-12) 
and Jeans (1919, sect. 18-23, pp. 20-6); the detailed demonstration of the necessity of an 
exchange of stability is found in in Poincare (1885). 

13 It might be possible to exclude them explicitly from the functional domain Q, avoiding 
any problems at the start. However, it is conceivable that such an exclusion would change 
the topology of Q and thus complicate the question of the connectivities of the space. For 
present purposes it is easier to deny them any place in the problem at the end. 
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a time-dependent, external tidal field, or any other case in which only part of an 
interacting system is modelled, the number of solutions cannot be determined from this 
development^. 



5. A Dynamical Example 



The simplest useful example of a dynamical system in astronomy is the two-body 
problem, dealing with a pair of bodies of reduced mass M in an orbit of total energy E 
and angular momentum J. Imposing a spherical coordinate system (r, 6, <p) with the orbit 
in the plane of the equator (8 = vr/2), the trajectory is given by 



1 + e cos ( 



(29) 



with e the eccentricity of the orbit and Rq = J 2 /GM. Defining the Jacobi functions in 
each of the coordinates as 5r — s, 89 — £, 5(p — r\ and the perturbations in energy and 
angular momentum as h and I respectively, one eventually finds 



dq 

d(p 

for e < 1, s 



£osin(0 

J r 
hGM 



(30) 
(31) 



2E 2 



F sin <b + e + 



El e 2 



e 2 + l 



Jh 



cos< 



e sin 



3e 2 



for e > 1, s 



2 1 + e cos i 
hGM 



sin 6 arctan 



2E 2 



F sin <b + e + 



El e 2 



Jh 



— tan — 
1 + e 2 



COS ( 



(32) 



e sm 



3e 2 



2 1 + ecos0 2y/e 2 - 1 



sin 6 In 



Ve 2 - ltan((/./2) + 1 + e 



Ve2^Itan(0/2) - 1 - e 



(33) 



where F is a constant used to adjust the zero point of s. The first expression for s is used 
for bound (elliptical) orbits, the second for unbound (hyperbolic). The practical difficulty 
of calculations using Jacobi functions is evident. 



14 This does not mean that Layzer's (1963) cosmic energy equation exempts all interesting 
distributions of astronomical objects from the results obtained here. An integral of energy 
still exists for any collection of masses interacting through gravity; Layzer's equation 
only states that a quantity based on comoving motions and coordinates, which resembles 
Newtonian energy in some respects, is not conserved. Since the number of solutions 
a problem has should not depend on which particular variables are used to write it 
down, results obtained herein using proper, inertial coordinates apply also to calculations 
performed in other ways. 
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The out-of-plane Jacobi function £ is, however, simple and gratifyingly general. For 
any eccentricity (indeed, even for unbound trajectories) conjugate points are found on 
diametrically opposite sides of the orbit. This is easy to picture: rotation of the orbit 
through an infinitesimal (or even larger) angle around a line from the orbiting body 
through the primary, certainly an allowed variation, leaves the opposite point unchanged. 

For very small e, that is for orbits close to circular, s becomes a simple sine function 
also, returning to zero after half an orbit. For e ~ 1, that is for orbits close to a parabola, 
analysis is a bit more complicated, though s is approximately sinusoidal and in no case 
does s reach zero again until after half an orbit. For e > 1, but not by much, s remains 
approximately sinusoidal. For very large e the approximation is better, that is the points 
where s vanishes are closer to being 180° in longitude away from each other; but the 
(unbound) trajectory might not include enough movement in longitude to provide a 
conjugate point for some initial points. 

The Jacobi function in longitude, rj, has a behavior which is in full even more 
complicated. However, note that its derivative is directly related to s. It can therefore 
not return to zero until well after s changes sign. Among the three Jacobi functions, 
then, £ has its first zero after exactly half an orbit, while the other two take longer; so 
The earliest zero for any perturbation in a two-body system occurs after half an orbit, so 
kinetic foci are 180° apart. 

Choquard's criterion is much easier to apply. There are two points where the 
momentum is normal to the gradient of the potential, at pericenter and apocenter; any 
pair of conjugate points must lie on opposite sides of one of these. Together with Morse's 
count of solutions to the Eulerian variational problem this means that there is an infinite 
number of solutions for a given set of end points, one for each half-integral number of 
revolutions of the orbit. As noted above, trajectories with energy lower than the lower 
of the potential energies of the end points are excluded from consideration. Those with 
positive energy have one least-action solution and possibly one saddle-point solution at 
finite times (depending on whether the first end point is taken far enough away from 
perihelion to allow the kinetic focus, 180° away in longitude, to appear on the trajectory); 
the rest at infinite times. 

Applied to systems with many bodies, saddle-point solutions correspond to some sort 
of multiple-pass trajectory. If there are two bodies in an orbit that approximates isolated 
two-body motion, they can generate kinetic foci for the whole system. 

Given a bound two-body system with a set of endpoints and a fixed time taken to 
go between them, the minimum-action solution will give a trajectory made up of less 
than half an orbit. The first stationary-action solution will contain more than half an 
orbit, a longer distance, which means a higher speed and thus higher kinetic energy. The 
second stationary solution will require at least three times the speed of the minimum 
solution, thus nine times the kinetic energy; few such orbits are bound. The situation 
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for a many-body system is rather more complicated, but for most astronomical systems 
a significant increase in the kinetic energy will make total energy positive and thus 
the system will become unbound. In this way the relatively small binding energy of 
astronomical systems severely limits the number of saddle-point solutions (unless there is, 
say, one or more tightly orbiting pairs of objects). 



6. Continuum Solutions 

The discrete body approach to galaxy dynamics is of course an approximation. It 
may be justified by the fact that present distances between galaxies are significantly 
larger than galaxy dimensions, or (more practically) on the basis of our ignorance of their 
detailed mass distributions (including such things as dark matter halos). But if we are 
to consider the motions of galaxies all the way back to their formation it becomes an 
increasingly bad approximation, and it would be better to consider a continuous fluid of 
gravitating matter. 

Indeed, the present picture of galaxy formation has them condensing out of a smooth 
fluid. It would be highly desirable to be able to follow this process in detail while requiring 
a certain configuration as a final end point. One could investigate, for example, the 
importance of mergers in galaxy dynamics, as well as the problems encountered by Dunn 
and Laflamme (1995) in matching a least-action calculation to an n-body simulation. 

However, in attempting this we are faced with a massive theoretical complication as 
the number of degrees of freedom goes from 3n to infinite^. Additional practical difficulty 
is involved with the increased complexity of the calculation, using three equations 
(continuity, Euler's and Poisson's) instead of one. However, it can be done, as Susperreggi 
and Binney (1994) have shown (though it tends to be computationally intensive). 

Consider, as a first approximation to a continuous-fluid situation, a large N-body 
calculation. Since the results of Morse Theory do not depend on the number of bodies, 
there still remains one minimum action solution and a finite number of stationary action 
solutions. (The bodies are now of all the same mass, and are labelled with, say, their 
ending coordinates instead of "M31"; but the Morse-based results are unchanged.) 
Adding more bodies increases the resolution of the simulation and the computational 
burden, but does nothing to the theory of solutions. Therefore, so far as a continuous 
fluid may be considered as made up of discrete masses, however tiny, there remains one 
least-action solution and one stationary solution for each possible value of the index. 



15 This is of less practical importance, as a continuum calculation always has some sort 
of short-wavelength cutoff (which is addressed in more detail below). 
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6.1. Orbit-Crossing and Kinetic Foci 



Giavalisco et al. (1993) identified orbit-crossing as a major cause of multiple 
solutions, that is, when trajectories from different parts of the fluid occupy the same 
point at the same time. This makes the mapping of velocity to distance (a major concern 
of observational cosmology) multiple-valued. However, our question — the number of ways 
the present velocity and density distribution can arise from the Big Bang — is different, 
and orbit-crossing is not necessarily relevant. 

To see this, consider a spherically symmetric part of a nearly uniform universe, of 
critical density for definiteness. Suppose that a small perturbation makes one shell slightly 
more dense than average and the shell contained immediately within it less dense. Over 
time the dense shell will expand at a slower rate than the universe as a whole, and the less 
dense shell faster; eventually their trajectories will meet, and there will be orbit- crossing 
(even with all shells expanding). 

To locate the kinetic foci, first write the dynamical equation of a shell which contains 
a mass m(r) within a radius r: 

r = -«^ (34) 



which has the Jacobi equation 



2 -^ps (35) 



r 

which, for shells near critical density, becomes 

s = Ars. (36) 
9t 2 v ; 

In any spherically symmetric case s can start from zero and go back to zero only after r 
changes sign. In a critical universe (and, indeed, in any universe before a Big Crunch) this 
never happens; thus there are no kinetic foci. Orbit- crossing does not necessarily generate 
kinetic foci. 

Now consider another nearly uniform universe, but this time allow several mass 
condensations to form. Place them in such a way as to generate two binary systems, and 
allow the tidal torque of each on the other to send them into bound orbits. In all this 
allow none of the trajectories of mass elements to cross. After half an orbit kinetic foci 
will be generated. Kinetic foci do not necessarily generate orbit- crossing. 

Certainly an orbit-crossing situation in the context of the cosmological problem 
demands that two mass elements start in the same place (where all mass elements start, 
the origin) and end in the same place (where their trajectories cross). At first glance this 
appears to involve two trajectories with identical (proper space) endpoints, and thus two 
solutions to the equations of motion. But a solution is made up of all the trajectories of 
the bodies included, and whether it is a saddle-point or a minimum is an attribute of the 
solution as a whole, not of any of these bodies. In fact the two bodies that share end 
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points in an orbit- crossing situation are two solutions to slightly different equations of 
motion, not two different solutions to the same equation. 



6.2. Potential Flow and Kinetic Foci 

A useful simplification, then, for a continuum least-action calculation would be one 
that eliminates closed orbits; that is, one in which there is no rotation. Susperreggi and 
Binney (1994) used a velocity field derived from a potential suggested by Herivel (1955): 

v(x,y,z,t) = Va(x,y,z,t). (37) 

The field thus derived is both laminar and irrotational; the first term refers to the fact 
that it can have no orbit-crossing, and the second to the fact that it can have no vorticity: 

V x v = (38) 

so they appear to have satisfied all parties. 

Unfortunately, it is possible to have rotation in a flow that has no vorticity. Equation 
(j38[) is satisfied by a velocity field whose longitudinal ((f)) component varies inversely with 
radius, oc R~ l ; Lynden-Bell has pointed this out and, moreover, shows that it is just 
the sort of field one expects from tidal interactions (Lynden-Bell 1996). A velocity field 
derived from a scalar potential can generate kinetic foci. 



6.3. Resolution and Kinetic Foci 

The number of solutions in a continuum calculation thus formally remains the same, 
even if the restriction to potential flow is imposed: one minimum and one or several 
stationary solutions. Considering the latter the situation can appear rather depressing, 
since any two-body orbit by any pair of mass-elements, no matter how small, will generate 
kinetic foci and thus multiple solutions. It seems somehow unfair that a cosmological 
simulation should lose its minimum status through half the orbit of its smallest binary 
star. In practical terms, this means that a continuum least-action algorithm which is 
strictly minimizing will find only one solution, the one without so much as a half-orbit, 
which is not necessarily the right one; while an algorithm which finds all stationary 
solutions will find many possible answers, with no clue as to which is more probable. 

But cosmological simulations rarely depict single stars. In practice there is always 
a scale below which no detail can be seen; kinetic foci on this scale cannot affect the 
minimum status of the calculation. In a very simple example, consider a triple star made 
up of one tight binary and one wide component. If all bodies are included, a solution will 
only be a minimum through half the orbital period of the close double. However if the 
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binary is modeled by a single mass, a solution will be a minimum through half the period 
of the wide component. It is a matter of choice which is the more important trajectory 
to calculate-or, alternatively, whether the computational burden of calculating several, 
perhaps many, stationary solutions is worth maintaining the higher resolution. 

In a more complicated situation setting the desirable resolution is also more 
complicated. In a rich galaxy cluster, for instance, the dynamical timescale of the center 
regions is much shorter than the outskirts, and varies continuously with radius. What 
particular scale is best for the calculation? The answer is not obvious. However, the 
question is not restricted to least action calculations, so it is at least a familiar one. 



7. Summary 

The important results of this study are as follows: 

// the action for the cosmological variational problem can be written in proper 
coordinates and an integral of energy exists, there is one minimum solution. Assuming 
Hamilton's Principle is used, there may be additional, stationary solutions, one for each 
number of kinetic foci, if multiple-pass trajectories exist. There is a finite number in 
total, limited by possible values of energy. Solutions containing at least one approximately 
two-body orbit which passes through more than 180° in longitude are not minima. 

Kinetic foci are reached only after the momentum is normal to the force for some 
body in the system. 

In so far as a continuous mass distribution may be approximated by an arbitrarily 
large number of individual masses, a continuum least-action calculation also has a single 
minimum solution, but generally a very large number of stationary solutions. These can 
be limited by setting a lower limit to the resolution of the calculation. The specific size of 
this resolution may be difficult to determine. 

A radial velocity, rather than a distance, can be used as an end point in a numerical 
variational calculation. Forms of the modified action required have been discovered by 
Giavalisco et al. (1993) and used by Schmoldt & Saha (1998). Using such an endpoint 
has no effect on the number or character of solutions. 

Orbit-crossing is not necessarily related to the number of solutions of a continuum 
calculation. 

It is a pleasure to thank Donald Lynden-Bell for drawing my attention to the 
least-action problem and suggesting the radial velocity action. I received valuable 
assistance in interpreting topological ideas from Wendelin Werner and Anthony Quas. 
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A. Variable End Points 



The following derivations follow Morse (1934) with some changes in notation and 
terminology. Courant and Hilbert (1953) have a derivation for the transversality condition 
which in fact results in the same formula; however, they require some assumptions about 
the end manifold which do not hold in the present situation. 



A.l. The Transversality Condition 



Suppose the problem to be that of minimizing the integral 

/= f 2 L{q h q h t)dt (Al) 
Jti 

subject to the condition that one or both of the end points are not fixed but must lie on 
end manifolds of some description. The solution to the problem is given by some extremal 
g = g(t). Admissible curves for the problem will be those with end points near those of g 
and which are continuous along with their first and second derivatives. These curves are 
described by the r functions ah(e) such that 0^(0) gives g. The end points in particular 
are given by 

t s = t s (ai, . . . , a r ) 

Qi = it ("l, • • • , Or) 
s E (1,2) 

(the superscript 1 or 2 refers to the initial or final end point). Observe 

qt(a h (e)) = qt(f( ah (e)),e) (A2) 
where h takes on the values 1 to r. 



Integral (|A1|) is now a function of e; considered this way, the first variation (by 
Liebnitz' Rule) is 



/'(e) 



dt s ' 



+ 



J i 



t2(e) y(^dq i+ 9Ldq i \ dt 

ti(e) Y \ d Qi de d( li 9e ) 



(A3) 
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After integration by parts and a bit of algebra, one obtains the Euler-Lagrange equations 
and 



. dL \ dt s dL dqi 

i 1 dq\ J de ^ dqi de 



1 2 



0. 



(A4) 



-i i 



Again, the parametrization by e is arbitrary. If de is multiplied out of the above equation, 
the normal form of the transversality condition is obtained. If the manifold on which 
the end point is allowed to vary is specified by means of the differentials dqf and dt s , 
(POD contains a condition fulfilled by the true minimizing end point. Conversely, the 
transversality condition can sometimes be used to gain some insight into the end manifold 
when only the Lagrangian and the fact of minimization are given. 



If the integral to be varied is changed from (|A1| ) to the velocity action, 
I* 



d I dL 

(lj dq) ]]l 



I - 



^ dL 



(A5) 



where j denotes those coordinates in which velocity rather than coordinate is fixed at 
the end point, the variation of the boundary term must be included in the transversality 
condition. A similar derivation to the above results in the velocity- action transversality 
condition: 



dL 
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A. 2. The Second Variation 



Applying Liebnitz' Rule again gives 
/"(e) = 



t2(e) d_ (dLdq^ dLdqA 
ti(e) de y 1 I % de dqi de ) 



+ 



+ 



(dL dqi dL dq\ \ dt s 
i \dqi de dqi de I de 



j i 



d dt s 

L/C- tXO 



+ 



d 2 t s 



de 2 



After some algebra this becomes 
/"(e) = 



ti(e) i \dedel Jti(e) ± de 2 \% dt \ dqi 



dt 



(A7) 
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+ 
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<9L d 2 gf 
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The second integral vanishes for extremals. 
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While this version of the second variation is useful, it may be made a manifestly 
symmetric quadratic form in the variations 

doth 



Uh 



de 



Using these in equation (AS) there results 
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For an extremal satisfying the transversality condition, the coefficients of d 2 ah/de 2 
as well as the last integral vanish; and we are left with the second variation integral, as 
in the case of fixed end points, and a symmetrical quadratic form in the variations at the 
end points. The variations at the end points and within the integral are related by 



dqj x , 

de ~ Y 
dqj 
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(All) 



where evaluation is carried out at the end points. Morse defines the quantities bhk for an 
extremal satisfying the transversality condition via 
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and uses this notation for his definitions of the index form. 
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If the velocity action is used, the second variation of the boundary term must be 
calculated and added to the expression above. Following the lines of the above derivation 
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one finds 
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Again, a symmetric form may be found for extremals which satisfy the transversality 
condition. Following the above derivation, one obtains 
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Fig. 1. — Poincare diagram connecting the two Least Action Principles. For a given 
set of end points (or manifolds) a slice of the functional domain, which by construction 
contains all the solutions to the problem, is plotted. The Hamiltonian action is the vertical 
coordinate, the Eulerian action the horizontal coordinate. For a given value of total energy 
(the Eulerian control parameter; say, E 1} E 2 or E 3 ) there will be a single minimum solution 
and a series of stationary solutions. For a given value of total time (the Hamiltonian 
control parameter, say t\, ti or £3) there will be at least one minimum solution. In terms of 
Catastrophe Theory, the minimum solutions form a chain of stable equilibria (shown as a 
solid line in the figure), the stationary solutions a chain of unstable equilibria (not shown in 
this figure for clarity). If there were two Hamiltonian solutions on any Eulerian branch of 
solutions, as shown here, it would require the meeting of two chains of Hamiltonian similar 
equilibria (at point C) without an exchange of stability. Such a situation is forbidden by 
Catastrophe Theory, as described in the text. 



